Skip to content

Conversation

@albertomercurio
Copy link
Contributor

@albertomercurio albertomercurio commented Oct 28, 2025

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

I add some specific methods for unary negation that may significantly simplify the expression when applicable.

Moreover, I have added a function that flattens the input Tuple of an AddedOperator. This avoids to have nested AddedOperators when just calling AddedOperator(ops).

If this pr will be merged, could a new version be released?

src/basic.jl Outdated
return ScaledOperator(L.λ, -L.L) # Try to propagate the rule
end
function Base.:-(L::MatrixOperator)
isconstant(L) && return MatrixOperator(-L.A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this allocates though

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, but this operation is usually performed once in the beginning when defining the operator. Then, the final operator is usually applied many times when solving an ODE or a linear problem.

This allows to have a more clean structure of the operator, which will be used repeatedly on a solver.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the ODE it will allocate for every step which takes a new Jacobian.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you make an explicit example? In my mind this is my workflow

some_f(a,u,p,t) = ...

# The allocations apply here, if the matrix is constant
A = ScalarOperator(rand(), update_function=some_f) * MatrixOperator(...) + ....

prob = ODEProblem(A, u0, tspan)
sol = solve(prob, Tsit5()) # It just updates the ScalarOperators. No allocations

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that's not using the WOperator

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You'd have to do a case where someone directly writes jac_prototype = MatrixOperator(J) I think, then it should allocate the extra terms instead of doing the lazy computation that is expected.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok now I'm doing

f = ODEFunction(brusselator_2d_loop; jac_prototype = MatrixOperator(float.(jac_sparsity)))
prob_ode_brusselator_2d_sparse = ODEProblem(f, u0, (0.0, 11.5), p)

function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = IncompleteLU.ilu(convert(AbstractMatrix, W), τ = 50.0)
    else
        Pl = Plprev
    end
    Pl, nothing
end

solve(prob_ode_brusselator_2d_sparse,
    KenCarp47(; linsolve = KrylovJL_GMRES(), precs = incompletelu,
    concrete_jac = true); save_everystep = false);

But still the same performances.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And memory usage and memory peak? The issue I'm worried about is that -A not being in-place will make there be a new temporary, and so it doubles the memory requirement on the system. I think this will oom larger cases.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm actually trying

Base.:-(L::AbstractSciMLOperator) = (println("HELLO"); ScaledOperator(-true, L))

without the other method definitions, just to check that it is actually calling it. But this is not. So maybe we should try another case?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok I have removed all the definitions that should allocate a new MatrixOperator, just to make it allocation-free.

@ChrisRackauckas ChrisRackauckas merged commit 73f73ac into SciML:master Nov 5, 2025
14 of 15 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants